Mergers of nonspinning black-hole binaries: Gravitational radiation characteristics 
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We present a detailed descriptive analysis of the gravitational radiation from black-hole binary 
mergers of nonspinning black holes, based on numerical simulations of systems varying from equal- 
mass to a 6:1 mass ratio. Our primary goal is to present relatively complete information about 
the waveforms, including all the leading multipolar components, to interested researchers. In our 
analysis, we pursue the simplest physical description of the dominant features in the radiation, pro- 
viding an interpretation of the waveforms in terms of an implicit rotating source. This interpretation 
applies uniformly to the full wave train, from inspiral through ringdown. We emphasize strong rela- 
tionships among the £ = m modes that persist through the full wave train. Exploring the structure 
of the waveforms in more detail, we conduct detailed analytic fitting of the late-time frequency evo- 
lution, identifying a key quantitative feature shared by the i = m modes among all mass ratios. We 
identify relationships, with a simple interpretation in terms of the implicit rotating source, among 
the evolution of frequency and amplitude, which hold for the late-time radiation. These detailed 
relationships provide sufficient information about the late-time radiation to yield a predictive model 
for the late-time waveforms, an alternative to the common practice of modeling by a sum of quasi- 
normal mode overtones. We demonstrate an application of this in a new effective-one-body-based 
analytic waveform model. 

PACS numbers: 04.25.Dm, 04.30.Db, 04.70.Bw, 04.80. Nn 95.30.Sf, 95.55.Ym 97.60.Lf 



I. INTRODUCTION 

The final merger of two black holes (BHs) having com- 
parable masses will produce an intense burst of grav- 
itational radiation, and is expected to be one of the 
strongest sources in the gravitational-wave sky. Mergers 
of stellar black holes are key targets for ground-based de- 
tectors such as LIGO, VIRGO, and GEO600, and knowl- 
edge of the merger waveforms is an important compo- 
nent of improving the dctectability of such systems. The 
space-based LISA detector will observe mergers of mas- 
sive black holes at high signal-to-noise ratios, allowing 
tests of general relativity in the strong- field, dynamical 
regime. 

Today, numerical relativity (NR) studies are beginning 
to progress toward a full description of black- hole binary 
merger systems. For noneccentric inspirals, this space is 
spanned by seven parameters: the symmetric mass-ratio 
77 = mim2 / {mi+m2)'^ , and the six combined components 
of the black holes' spin vectors. Considerable study has 
been focused on the fiducial center point of this parame- 
ter space, the case of equal-mass nonspinning black-hole 
mergers. After the series of breakthroughs that ushered 
in an era of rapid progress in the field H, d, [1] j sev- 
eral investigations assessing the accuracy of the available 
equal mass waveforms and applying them to data analy- 
sis were conducted [1, H, 0, E, iSl] ■ 

In this paper, we undertake a descriptive study of the 
waveforms generated in the late inspiral and merger of 
black-hole binaries for the subspace of nonspinning black 
holes, parametrized only by 77. Our study is based on a 
series of numerical simulations, discussed in Sec. lIIII cov- 



ering at least the last ^ 4 orbits of nonspinning black- 
hole binary mergers with mass ratios extending to 6:1 
(?7 w 0.122). Several of the simulations presented here 
have already been applied in a recent paper, focusing on 
the development of a faithful analytic waveform model 
p^ . Here we provide details of these and additional 
simulations, together with considerable analysis, focused 
on providing a qualitative and quantitative picture of 
how the waveforms from nonspinning black-hole merg- 
ers depend on 77. Nonspinning black-hole binar y m erger 
waveforms were previously examined in Ref. [ll|, but 
our analysis is novel and complementary to that work. 
Our descriptive presentation puts emphasis on the re- 
lationships between waveforms from the different mass- 
ratio cases and different harmonic modes, with references 
to Ref. [llj where related observations have been made. 
Our approach to describing the inspiral-merger-ringdown 
transition is particularly distinct, founded in a uniform 
approach that describes all stages of this process in sim- 
ilar terms, and ultimately suggesting a complementary 
physical picture. 

Black-hole-binary merger waveforms have been noted 
for their "simplicity." For the nonspinning systems the 
simple physics of the coalescence is exposed by a spheri- 
cal harmonic decomposition of the waveforms. In Sec. lIVI 
we walk readers through the basic features of the radi- 
ation, characterizing amplitude and phase evolution of 
the multipolar components, and discussing relationships 
among the simulations representing different mass ratios, 
and among the multipolar components of each simula- 
tion. As we analyze the waveforms we develop a concep- 
tual interpretation of the basic waveform features. In this 
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interpretation we consider the structure of an implicit ro- 
tating source, which could have generated the measured 
radiation through its rotational motion. This allows a 
uniform interpretation that applies throughout the coa- 
lescence process: inspiral, merger and ringdown. 

In Sec. |Vl we examine the strong final burst of radi- 
ation beginning ^ 2QM before the formation of a com- 
mon horizon. We quantitatively describe the phasing in 
terms of an analytic model, based on a continuous, mono- 
tonically increasing frequency. We find, in particular, 
that the peak rate of change in frequency, appropriately 
scaled, is the same across all £ ^ m modes and mass 
ratios. We also identify relationships among the mode 
amplitudes and phases, which are connected to an ap- 
proximately linear relationship between angular momen- 
tum and frequency: cPj/duP' ^ 0. We interpret these 
relationships in terms of the implicit source. 

Finally, In Sec. IVIl we demonstrate the utility of what 
we have learned in our waveform characterization by ap- 
plying some of the quantitative features we have uncov- 
ered in a new variation on the analytic waveform model in 
[l0| . which was based on the Effective-One-Body (EOB) 
resummation of the Post-Newtonian(PN) approximation 
to inspiral dynamics [l^ . In particular, wc provide a dis- 
tinct late-time waveform model, alternative to the com- 
mon "spectroscopic" model [ij, [I3| based on sums of 
quasinormal mode overtones. 

II. OVERVIEW 

We begin with some examples of gravitational strain 
waveforms as they might be observed by gravitational- 
wave instruments. In observational work, and PN analy- 
sis, it is customary to describe the radiation in terms of 
gravitational-wave strain, h. In representing the strain, 
it is convenient to combine the two real waveform polar- 
ization components, hj^ and /ix, hito one complex strain 
waveform, 

/i = /i+ + i/ix ■ (1) 

We decompose the strain waveforms measured on a 
sphere of radius i?oxt, into spin- weighted spherical har- 
monic components, him- The details of the decomposi- 
tion, and how the waveform information is extracted from 
the numerical simulations, are given in Appendix \K\ 

The waveforms in this section are aligned in time and 
phase so that the maximum strain amplitude occurs at 
t = 0. The remaining figures of this paper will be aligned 
in a similar way, but with t = Q marking the time of peak 
(2.2) mode energy flux, E22 (unless stated otherwise). 

Fig. [T] shows waveforms from mergers of nonspinning 
black holes for various mass ratios, as observed at dis- 
tance R on the rotational/orbital axis of the system. 
The figure shows h+ for each of the four mass ratios 1:1, 
2:1, 4:1, and 6:1. For these observers the observed wave- 
forms will be circularly polarized, so that /ix is 90° out 
of phase with . We use units in which G = 1 and c = 1 



and express both time and spatial distances in terms of 
the total mass Af, where Af ~ 5 x 1O~^(M/M0)sec ~ 
1.5(M/M0)km. 

More typically, the observer will not be located on the 
system's orbital axis. The left panel of Fig. [5] shows 
/i+ for the 4:1 case. The strain is measured at an az- 
imuthal angle of 0° and various inclinations^. The de- 
tailed shapes of the waveforms change as the system is 
reoriented so that the observer moves off the system's 
rotational axis. For larger inclinations (closer to being 
viewed edge-on) there are notable modulations at half 
the base gravitational-wave frequency. 

The right panel of Fig. [2] shows for different mass 
ratios oriented at an inclination of 90° and an azimuthal 
angle of 0°. For this orientation, constitutes the full 
strain waveform. For larger mass ratios, the lower fre- 
quency modulation increases. For gravitational-wave ob- 
servations of sufficiently strong binary black-hole sources, 
the types of differences shown in Fig.[2]could be exploited 
to estimate the inclination and mass-ratio of the source 
system. 

For observational purposes, the combined waveform in- 
formation encoding all possible source orientations can be 
conveniently represented in terms of spin- weighted spher- 
ical harmonic components [see Eas. (jA2p - (|A4p ]. providing 
a neat description of the leading waveform features. The 
multipolar decomposition is even more valuable as a tool 
for exposing the hallmark simplicity of the merger ra- 
diation. The readily apparent simplicity in the waves 
viewed from the system's orbital axis in Fig. [T] extends 
to each of the spherical harmonic components. Viewed 
off axis, these components linearly combine to yield the 
more complex appearance of the waveforms in Fig. [2] 

This characterization of the gravitational radiation 
from a merging black-hole binary in terms of circular 
polarization was first recognized in the Lazarus project 
studies [IBl- In this picture, the radiation can be repre- 
sented by a slowly varying amplitude and a polarization 
angle; sec Eq. ^ below. This description relies on how 
the radiation appears to distant observers located on the 
rotational axis of the system. Other observers will typ- 
ically see elliptically polarized waves, having a generally 
simple pattern that conforms to the rotational nature of 
the source. In the equatorial plane, the radiation reduces 
to the plus polarization, corresponding to the observer 
seeing no circulation in the source. Looking along the 
negative rotation axis, the observer sees circular polar- 
ization with the opposite hclicity. 

Each of the spherical harmonic waveform components 



The inclination angle is defined here as the angle between the 
linc-of-sight with respect to the detector and the orbital axis of 
the binary. This is the same angle referred to as "inclination" in 
the PN/NR literature, and most equations are constructed using 
that definition. However, the astronomical literature has often 
defined inclination to be the angle between the line of sight and 
the orbital plane of the binary, resulting in a 90° inconsistency. 
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FIG. 1: The plus component of the strain, h+, is shown for different mass ratios using the £ = 2 and £ = 3 modes. The observer 
is located at distance R along the axis of the system, at inclination (footnote [TJ 9 = 0° and azimuthal angle = 0°. The strains 
are scaled by symmetric mass ratio and aligned such that, for each mass ratio, the peak of /i22 occurs at f = 0. The phases are 
rotated such that the phases are 0° at t = 0. 
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FIG. 2: Waveform variation with inclination. The left panel shows how the wave shape changes with inclination 6 for a 4:1 
mass ratio. h+ is plotted, using the sum of multipolar modes up to ^ = 5. The curves are labeled by the inclination in degrees 
of the observer, and the time axis is labeled so that the peak of /i22 occurs at t = 0. The right panel shows how the wave 
shape changes with mass ratio. The waveforms h+ are computed at an inclination of 90° from the £ = 2 and £ = 3 modes. The 
different mass ratios are aligned in time so that the peak of /122 occurs at t = 0. The rotational phases are matched to 0° at 
t = 0. In both panels, the azimuthal angle is 0°. 



exhibits circular polarization with steadily varying phase 
and amplitude, providing a natural framework for devel- 
oping a practical and intuitive understanding of binary 
black-hole merger radiation. Our basic waveform descrip- 
tion, in Sec. IIV| and the more detailed analysis that fol- 
lows, is based on this spherical harmonic decomposition. 

As wc describe the waveforms, we will also suggest 
a heuristic interpretation of what the radiation tells us 



about the motion and structure of the binary black-hole 
source. In the weak-field description of radiation from 
a rotating object, the multipolar waveform components 
of the gravitational radiation can be associated with dy- 
namics of multipolar moments of the radiating source 
It is useful, in conceptualizing the full coalescence 
radiation from inspiral through merger and ringdown, to 
think of the multipolar radiation description as provid- 
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ing information about the motion of a changing source 
object, described as a sum of several multipolar mass 
moments. This source object is what we will interpret as 
an effective rigid rotator radiation source, with a slowly 
changing structure. We refer to this as the implicit ro- 
tating source (IRS). 

In the process of coalescence, the source begins as a 
separated black-hole binary system and ends as a single 
distorted black hole. For nonspinning binary mergers, 
numerical and PN results consistently indicate that the 
radiation is circularly polarized, in the sense first rec- 
ognized in the Lazarus project studies (l5| . not only in 
the inspiral, but uniformly through the merger and ring- 
down. In our conceptual source description, this pattern 
of circular polarization is consistent with radiation gen- 
erated by rotational motion of each source multipolar 
moment, where the polarization phase is tied to the in- 
stantaneous orientation of the source. Similarly, we think 
of the amplitude of the radiation multipole as related, 
through some generalization of the quadrupole formula, 
to the magnitude of the implicit source multipole. 

We can write each multipolar component in a specific 
polar form natural for circularly polarized radiation: 



rhi„i (t) 



im*<|2 (*) 



(m > 0) ; 



(m<0). 



Each amplitude Him is expected to be a slowly varying 
function of time and can be conceptually considered as 
a function only of the magnitude of the source multipo- 
lar moment and its rotational frequency. The additional 
sign for some m < cases allows a consistent interpreta- 
tion for the component phases and the component am- 
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plitudes. The waveform phase is given here by to$ 

but in most of our analysis we refer to which we 

call the rotational phase. In our terms of our implicit 
source heuristic, the rotational phase for some particular 
(£, to) mode can be thought of as the azimuthal orienta- 
tion of the specified multipolar source component. In the 
inspiral, where the source can be considered as a sepa- 
rated binary, <i>^^' should coincide with the orbital phase, 
independent of (£, to). In the post-Newtonian expansion 

all the 'I'^^ defined here agree to at least 2PN order, 
while the amplitudes Him remain real and non-negative. 
The equatorial-plane symmetry of these mergers ensures 



that Hf. 



lira = i?^(-m) and that - 

need consider only m > modes for this analysis. 

The expansion Q is not appropriate for the m = 
modes. This points to an important caveat to our implicit 
rotating source interpretation of the radiation, that it ap- 
plies to the degree that the radiation is circularly polar- 
ized. While not strictly vanishing, the m = waveform 
components, and other deviations from circular polariza- 
tion are generally extremely small, and largely unmea- 
surable at the resolutions of the simulations we study. 
For the most part, we will not address deviations from 
circular polarization, and the consequent limitations of 



(h) 



so that we 



our implicit-rotating-source interpretation in this paper, 
focusing for now on the dominant features of the radia- 
tion. 



III. SIMULATIONS 

Our analysis is based on four simulations, representing 
mass ratios 1:1, 2:1, 4:1 and 6:1. Results from the 1:1 and 
4:1 simulations have appeared in previous publications 
(01 and [l3|, respectively). More recently, higher accu- 
racy simulations have been presented by other groups for 
the 1:1 case [1, [l^. Our older waveform is sufficiently 
accurate for our present purpose, to examine the general 
features of the nonspinning merger waveforms. 

Our numerical simulations are carried out with 
the HAHNDOL evolution code [l^, which uses finite- 
differencing methods to solve a 3-1-1 formulation of Ein- 
stein's equations on a Cartesian grid. For initial data 
we solve the elliptic equation given by Brandt and 
Briigmann for conformally flat data in which the black 
holes are represented by punctures [1^. This is per- 
formed numerically using the multigrid solver AMRMG 
[2l[, which is second-order-accurate but tuned to give 
truncation errors typically much smaller than those pro- 
duced by the evolution code. The momentum param- 
eters are chosen according to the 2PN-accurate quasi- 
circular approximation given by Kidder p^ . which has 
been found to result in low eccentricity. We evolve these 
data using the moving puncture method [1, [2^ with 
a modified version of the Baumgarte-Shapiro-Shibata- 
Nakamura equations [l^, . Specifically, as suggested 
in [2^ , we replaced the conformal factor variable (j) with 
exp(— 20), which vanishes at the punctures. Further, we 
added the constraint-damping terms suggested in [23], 
and the dissipation terms suggested in j28l . [29j . For 
the gauge we use the specific l4-log lapse and Gamma- 
freezing shift conditions recommended for moving punc- 
tures in (soj . 

Accurate simulations require adequate spatial resolu- 
tion near the black holes (length scales ~ M) as well 
as in the wave zone where the gravitational waves are 
extracted (length scales (10 - 100)M). To this end, 
the grid has multiple refinement levels, determined adap- 
tively near the black holes, but fixed in regions farther 
away (typically, |a;| > 30Af) where the waves are ex- 
tracted; all grid refinement is handled within the frame- 
work of the software package paramesh [31| . The adap- 
tive mesh refinement criterion near the black holes is 
designed to keep the scale of the square root of an in- 
variantly defined curvature component, known as the 
Coulomb scalar [s^ . [ssj . roughly constant with respect 
to the grid spacing. Interpolation in guard-cells between 
refinement regions is fifth-order-accurate, coupling with 
differencing stencils to yield at least fourth-order accu- 
racy in the bulk. 

Spatial derivatives are taken by sixth-order-accurate 
differencing stencils, with the exception of advection 
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derivatives, which are handled by fifth-order-accurate 
mesh-adapted differencing for greater stabihty [33|^. 
Time integration is performed with a fourth-order Runge- 
Kutta algorithm. 

The initial configurations of the simulations we analyze 
are given in Table [H In each case, the initial separation 
was chosen to be large enough to result in at least five 
orbits. The finest resolution, hf, ranged from M/32 to 
3M/224, as required to adequately resolve the black hole 
with the smaller mass in each case. The outer boundary 
was typically at |a;| > lOOOM, far enough away to pre- 
vent refiections from reaching the wave-extraction region 
during the simulation. 



We have measured the individual black-hole masses 
using the apparent-horizon mass mi, a quantity calcu- 
lated from the area of each hole's horizon, which we lo- 
cate using the AHFinderDirect code HH]. From these 
horizon masses, we calculate the symmetric mass ratio 
r] = mim2/{mi + 7712)^- This gives the most precise 
specification of the actual mass ratio attained in our sim- 
ulations. In the text we will refer to the simulations by 
the mass ratio (e.g. 4:1). 

We define the total, infinite-separation, mass M of the 
system as an analogue for the total rest mass parameter 
used in PN studies. We measure M in two ways: 

MAH = mi+m2, (3) 

the sum of the individual BH horizon masses, and 

A/rest = A/adm - Eb, (4) 

defined as the difference between Madm , the total energy 
of the initial data, and the (negative) binding energy of 
the binary. The binding energy Eh is estimated from an 
effective-one-body PN treatment [l^, given the initial 
angular momentum Jq. The result shows a very close 
correspondence between Mah and Mrest j with differences 
at the level 10"''. For the rest of this paper, we use 
M = Mah, except for the 1:1 simulation data, where 
A/ah was not available for technical reasons. 

In interpreting the late-time radiation, it is valuable 
to know the mass and spin of the final Kerr black hole 
formed by the merger. We discuss the state of the final 
black hole, determined consistently by several means, in 
Appendix [Cl 

For the 4:1 mass ratio case, we have carried out runs at 
three different resolutions in order to assess the quality 
of the simulations. The convergence of the constraints 
and waveforms is discussed in Appendix |B] 



^ Sixth-order center-difTerencod advection is unstable, and sixth- 
order lopsided advection is too costly in terms of PARAMESH 
guardcoUs, which motivated our particular modification. 



The most important products of our simulations are 
the gravitational radiation waveforms, which we extract 
from the evolved simulation data as explained in Ap- 
pendix [Xj Strain-rate waveforms for the 4:1 case at vari- 
ous resolutions are shown in Fig. [3l where the times and 
phases have been shifted to agree at the moment of peak 
energy flux, as is generally done in our analysis below. 
We can get some measure of the error in the waveforms by 
comparing the difference between the high and medium 
resolution simulations. Fig. [3] shows the relative differ- 
ences in amplitudes, scaled by the high-resolution result. 
Ignoring the high frequency noise, the (2,2)-mode differ- 
ences (upper panel) indicate a combination of a secular 
amplitude difference and a sinusoidal effect, which results 
from the combination of the eccentricity in the orbital 
dynamics and the difference in peak time due to limited 
resolution. These combine to give differences generally at 
the 3% level, somewhat smaller at late times. The eccen- 
tricity plays less of a role in the (4,4) differences (lower 
panel), as the relative secular error is much (~ 5 times) 
larger. Sinusoidal eccentricity effects are also visible in 
the phasing error [see Fig.[S]. Overall, we find waveform 
amplitude and phase errors to be consistent with between 
fourth- and fifth-order convergence [see Appendix IB] . 

Assuming fourth-order convergence, and using 
Richardson extrapolation, our nominal expectation for 
these simulations leads to an error estimate for the high- 
resolution simulation applied in our analysis of ~ 1.2 x 
the difference shown in Fig. 2) To be conservative, we 
could instead assume second-order convergence, which 
would lead to an error estimate of ~ 2.8 x the difference 
shown in Fig. ID 

The errors for the 2:1 case should be comparable to the 
4:1 case. The resolution for the 6:1, scaled by the smaller 
black hole's mass is about 15% lower than lowest resolu- 
tion 4:1 simulation, suggesting errors eight times larger, 
if we conservatively assume fourth-order convergence and 
that the errors around the smaller black hole dominate. 
The errors for the 1:1 mass ratio case are discussed in 
Ref. 0. 



IV. DESCRIPTIVE ANALYSIS OF 
WAVEFORMS 

In this section we provide a descriptive analysis of the 
waveforms from our simulations. We try to serve two 
purposes in analyzing the radiation. In the first place, 
we are hoping to provide material for gravitational-wave 
observers, and others outside the field of numerical rel- 
ativity, which makes clear some of the general charac- 
teristics of the radiation from these mergers. Beyond 
that, we also push the analysis in more detail, hoping 
to generate deeper insight into the physics which gener- 
ates the radiation. Through this analysis we explore the 
similarities and differences for the various mass-ratio sim- 
ulations, and among the different multipole components 
of the 4:1 case case. In this way, we examine the wave- 
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TABLE I: Physical and numerical parameters of the initial data for all the runs presented. mi,p and m2.p are the puncture masses 
of the two pre-merger holes, ro and Po are the initial coordinate separation and (transverse) linear momentum, respectively, 
giving rise to a total initial orbital angular momentum Jq. hf is the spatial resolution of the highest refinement level for each 
run. Madm is the total energy of the initial data. The total, infinite-separation mass M of the system is measured in two 
ways " Mah, the sum of the initial (apparent) horizon masses of the two holes, and Mrest, the sum of the ADM energy and the 
binding energy from efi^ective-one-body theory [T^]. Finally, 77 is the resulting symmetric mass ratio, as determined from the 
two holes' horizon masses. 



Mass ratio hf 


mi^p m2,p ro Po Jo 


Madm Mah Mrost r] 


1:1 M/32 


0.4872 0.4872 10.800 0.09118 0.9847 


0.9907 ■ ■ ■ 1.0005 0.2500 


2:1 3M/160 


0.3202 0.6504 8.865 0.09330 0.8271 


0.9889 0.9989 0.9990 0.2228 


4:1 3M/224 
Af/64 
3A//160 


0.1890 0.7900 8.470 0.06957 0.5893 
0.1890 0.7900 8.470 0.06957 0.5893 
0.1890 0.7900 8.470 0.06957 0.5893 


0.9929 1.0003 1.0004 0.1601 
0.9929 1.0003 1.0004 0.1601 
0.9930 1.0003 1.0005 0.1601 


6:1 M/64 


0.1338 0.8490 8.003 0.05559 0.4449 


0.9942 1.0000 1.0001 0.1226 




-100 
t-t ,(M) 

peak(l,m) ^ ' 



FIG. 3: The (2,2) (top panel) and (4,4) (bottom panel) strain- 
rate waveform modes from the 4:1 mass-ratio case. Three 
resolutions are shown, time- and phase-shifted to match at 
their respective peak amplitudes. Excellent agreement is seen 
to persist throughout most of the simulation. 
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FIG. 4: The relative amplitude 6rrOr (-Amcdium ~ nno 
for the (2,2) (top panel) and (4,4) (bottom panel) strain-rate 
waveform from the 4:1 mass-ratio case. 



study the waveform amplitudes and the associated en- 
ergetics of the merger in IIVBI In IIV CI we address the 
polarization phase of several waveform modes, relating 
them to a common implicit source phase. Next, in SccfVl 
we will examine the late-time frequency evolution and the 
relation to amplitude in more quantitative detail through 
the inspiral-merger-ringdown transition. 



Strain-rate 



form amplitudes and energy, and the waveform phasing. 
As we proceed, we will interpret the results in terms of 
our implicit-rotating-source model, building up a heuris- 
tic description that applies through the inspiral, merger 
and ringdown of the binary. 

Following a brief discussion of strain rate in llV Al we 



In the Introduction, we motivated the spherical har- 
monic phase and amplitude waveform decomposition 
with a discussion in terms of strain h. In analyzing our 
numerical simulation results, however, we can work more 
directly with the strain rate h(t) = dh(t)/dt [see Ap- 
pendix 1^ for a more detailed discussion] . As with the 
strain decomposition we will expand the strain rate 
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FIG. 5: The absolute phase error (v'mcdium — V'flnc) (in radians) 
for the (2,2) (top panel) and (4,4) (bottom panel) strain-rate 
waveform from the 4:1 mass-ratio case. 
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FIG. 6: The (2,2) mode of the peak power versus mass-ratio 
rj (solid circles), with a numerical fit to the quadratic-quartic 
formula (|6]). Also shown is the full peak power, divided by 
two (asterisks). The growing difference between the two as 
77 ^ indicates the increasing importance of the subdominant 
modes in unequal-mass cases. The right-hand scale gives the 
peak power in c.g.s units. 



(m > 0) 



-iA£™e™*^"« (-1)^ (m<0), 



(5) 



with Aim real and non-negative.'^ Direct differentia- 
tion of (21) reveals the relationships between the phases 
and amplitudes defined in ([2]) and ([5]), with Aim ^ 



Um + 0{5PN) while differs from only 
at 2.5PN order. Note that the differentiation produces a 
phase shift of tt/2, so that strain-rate phases should be 
defined in the slightly unconventional form ^ if we wish 
to preserve the property that all are equal to orbital 
phase in the limit of well-separated binaries. This rep- 
resentation allows a more meaningful comparison of the 
phases of different multipolar modes, ensuring that in- 
stantaneous phase corresponds to the orientation of the 
binary system during the inspiral. 

The strain-rate amplitude is directly related to the ra- 
diation power for the {£,m) mode by Eq. (jATp . ~ 
(ylfm)^/167r. Henceforth we shall use the strain-rate- 
derived rotational phase ^im, rather than We also 
limit our presentation to the 7ti > modes, as equatorial 
symmetry implies ^ern = $f(-m) and Ai,n = 

Unless otherwise indicated, in the remainder of the pa- 
per, the time axis of each plot will be shifted so that the 
peak of £^22 (and hence of the strain-rate amplitude A22) 
occurs at t = 0. As the (2,2) mode is strongly dominant, 
this will closely approximate the peak time of the total 
E. 



B. Amplitude and energetics 

We first study wave amplitudes across modes and mass 
ratios. Since we are examining strain rate waveforms, the 
modal energy flux is effectively equivalent to the square of 
the mode amplitude, as in Eq. (jA7p . Preferring the most 
physical language we will express our modal amplitude 
comparisons as energy flux comparisons. In terms of the 
implicit rotating source model, we can think of the energy 
carried by the radiation as energy lost by the source. 

In Fig. [6l we plot the actual peak values of the dom- 
inant (2,2) energy flux contribution from Eq. (|A7P as 
a function of symmetric mass ratio rj. As this mode 
contribution is proportional to |/i22p, and we expect 
lim^^o \h22\ — * (the test-particle limit), we fit it to 
a quadratic-quartic form, obtaining the fit^: 



E22iv) = (4.40 ± 0.17) X 10- V 



-(5.43 ±0.31) x IQ-^T]^. 



(6) 



We also plot the peak of the total energy flux for each 
mass ratio, scaled by one-half, since the (2,2) and (2,-2) 
modes contribute equally to E. The difference between 
^'TOTAl/2 and E22 increases as 77 — > 0, reflecting the 
increased importance of other modes for unequal masses. 

Aside from the value of the radiation maxima, it is in- 
teresting to see how the radiation power evolves in time 
near the peak. In Fig. [3 we show shapes of the domi- 
nant (2,2) contributions to the peak energy fluxes (|A7[) 



^ A similar expression, differing from ((2J by an overall sign, would 
be equally appropriate for direct interpretation of numerically 
derived 1/14 waveforms. 



* A quadratic-cubic form is equally plausible, but fit the numerical 
data worse in this case. 
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for each mass ratio. These are scaled to the same peak 
height to allow shape comparison, and shifted in time so 
that the peaks of E22 arc aligned. We note the strik- 
ing similarity of the peak shape and duration across all 
mass ratios. Dm^ing the late-inspiral phase, the more ex- 
treme mass ratios appear to radiate more energy; how- 
ever, since we have normalized each curve by peak height, 
this only means that the equal-mass binary experiences 
a steeper climb to its peak power rate. Nevertheless, the 
different mass ratios follow similarly shaped tracks ap- 
proaching merger, and differences have been effaced by 
~ lOM before the peak power. The post-peak portion 
of the curve is determined by the dominant quasinormal 
mode (QNM) damping time, which varies only slightly 
with underlying Kerr spin for moderate spins [see, for in- 
stance, tables in The lower spins on the black holes 
formed by smaller-77 mergers should cause power to fall 
off faster in these cases. The inset in Fig. [7] shows the 
difference in fall-off rate relative to the equal-mass case 
from 20M after the peak. 

In Fig. m we concentrate on the 4:1 mass ratio, plot- 
ting the strongest modal contributions to the fLux, again 
scaled to the same peak height. Here we have applied 
only an overall time shift, so that the total flux peak is 
at t = 0. The energy profiles of the radiation burst in 
modes with £ = m are similar, all peaking at approxi- 
mately the same time. The subdominant modes are rel- 
atively stronger in the burst than in the inspiral, so that 
they show up here as weaker in the approach to the peak 
even after scaling to the peak radiation. We note that 
the subdominant modes with £ = m, (3, 3) and (4, 4), are 
particularly similar in this regard. Similarity in £ = m 
modes, and distinction in the other modes is a general 
feature of the bursts in several ways. 

The shape of the peaks with i m are particularly 
distinct. The (2,1) mode peaks particularly late, and the 
burst is much stronger than the inspiral. We note that 
the (3,2) mode shows a double-bump in its contribution 
to the energy flux. From Fig. O this appears to be robust 
in its gross shape over resolution and extraction radius. 
From Fig. [TO]we note, however, that the extent of this 
double-bump effect is very dependent on mass ratio; it 
is less in the 6:1 case, and not evident at all in the 2:1 
case. Later we will also note irregularities in the late-time 
frequency evolution of this mode, apparently indicating 
a deviation from circular polarization in this case. 

In both Fig. [7] and Fig. [5] wc have normalized the 
modc-flux peaks for the purposes of shape comparison. 
It is also important to understand the relative strengths 
of each mode. Wc show in Fig. [11] the relative mode con- 
tributions to E(t) for several dominant modes over the 
final inspiral and merger of the 4:1 case. In the merger- 
ringdown peaks, as in the inspirals, the £ = m modes 
dominate the energy flux, followed by the £ = |m| -I- 1 
modes. More discussion of the relative mode strengths 
for general nonspinning mergers is given in Ref . (ll| . 

In addition to the direct mode contributions, wc plot 
the PN-derivcd energy "partitions" - the power emitted 




t(M) 

FIG. 7: The shape of the (2,2) contribution to E{t) for all 
mass ratios, scaled to unity at the peaks. We have shifted the 
time axis so that the peak in £(2,2) occurs at f = 0. The inset 
focuses on the ringdown portion of the curves. 




t(M) 

FIG. 8: The shape of the different multipolar contributions 
to E(t) for the 4:1 case, scaled to unity at the peaks. We 
have shifted the time axis so that the peak in £{2,2) occurs at 

in each mode as a fraction of the total, according to the 
leading-order "restricted" FN expressions found in Eqs. 
(30)-(36) of [l3|, where the underlying orbital frequency 
was derived from the (2,2) mode. These are shown in Fig. 
[TT]using dashed lines. Until near the peak time, where we 
discontinue the FN curves, we see that the partitioning 
tracks the numerics very well for all modes except (2,1), 
which grows visibly faster than the restricted amplitude 
prediction after ~ lOOM before merger. A similar study 
jl(j | showed that the restricted (leading-order) approxi- 
mation for the amplitudes consistently overestimates the 
strength of the radiation. This shows that the leading- 
order partitioning of energy can provide a simple, but 
more accurate, approximation of the mode amplitudes. 
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FIG. 9; The (3,2)-mode "double- bump" in E{t) for the 4:1 
case. We show the inner two extraction radii {2QM and 45Af) 
for each of our three resolutions, shifting the time axis so that 
the peak in E(^2,2) occurs at f = 0. Though there is variation 
in the detailed shape, the overall double-bump envelope seems 
to be robust. 
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FIG. 10: The (3,2)-mode in E{t) for the 2:1 (red/dashed), 4:1 
(green/dot-dashed), and 6:1 (blue/dotted) cases, extracted at 
-Rcxt = 45M. We have shifted the time axis of each data set 
so that the peak in -E(2,2) occurs at f = 0. The amplitudes are 
unsealed. It is clear that the 4:1 case has the most pronounced 
deviation from a single well-defined peak. 



We will take advantage of this in Sec. lVII to provide more 
accurate amplitudes in a variation on the analytic EOB- 
based waveform model studied in (lo| . 



C. Waveform phasing 

In gravitational-wave observations, the waveform 
phase provides most of the time variation in the signals, 
and consequently is critically important in encoding ob- 
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FIG. 11: Relative power -Ef,m for different multipolar modes 
of the 4:1 case, on a logarithmic scale (solid lines). For com- 
parison, we also plot (dashed lines) the fractional power ex- 
pected from each of these modes based on E g. (IA6I) . and the 
restricted PN waveforms of Eqs. (30)-(36) of [l3~The differ- 
ent modes have been aligned to peak at Z; = 0. The right-hand 
scale gives the power in c.g.s units. 



servablc information about the source. Here we consider 
the phasing of the leading spherical harmonic waveform 
components. Following the discussion above, we concep- 
tually interpret each waveform phase as describing the 
orientation of a particular (^, m) multipole of an implicit 
rotating source of the gravitational waves. 

Direct comparative analysis of phases provides a 
stronger probe of the phase relations among the multi- 
polar modes than the comparative analysis of frequencies 
conducted in a number of previous studies of numerical 
simulations. Here we will discuss waveform phasing in 
terms of as it appears in ([5]), which we will compute 
from each strain-rate mode him(t). As noted in Sec. [Ill 
we expect all phases to agree in the large-separation limit. 

We first compute the strain-rate waveform phase ipim 
using the conventional decomposition: 



(7) 



with Vtrn real and non-negative. Then, setting this equal 
to (O, and solving for the rotational phase ^em{t), we 
find 



^ {'^(rnit) - f + ^Tniim) [m > 0) 

^ {fem{t) + ^ + 2nnem + (tt) (m < 0) , 



(8) 

The ^ term results from the factor i in ([5|), while the 
(—1)^ factor there produces the in term for to < 0. The 
nirn terms express the 2tt ambiguity in defining wave- 
form phase. Considering any {£, m) mode in isolation 
leads to an TO-fold degeneracy in the associated rota- 
tional phase. We resolve this degeneracy by choosing the 
pair {n22,?^33} G {(0,1) x (0,1,2)} that yields the clos- 
est consistency between $22 and $33 at early times (near 
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t = — 400M). We then determine the remaining n^m for 
closest early consistency with $22- This gives us a phase 
for each mode that can be interpreted as the rotational 
phase of the implicit rotating source that produced that 
component of the radiation. 

Fig. [12] shows the rotational phase ^im from several 
modes of the strain-rate waveforms from the highest- 
resolution (3M/224) 4:1 run, together with the rotational 
phase calculated from the tracks of the punctures. The 
left panel shows that for the inspiral portion of the evo- 
lution, all waveform phases, extracted at i?oxt = 45M, 
agree extremely well, except for the (2, 1) and (3, 2) 
modes, which differ by a significant part of a radian. The 
relative difference of each mode from the (2,2) mode is 
shown in the right panel. The differences between the 
£ = m modes are ^ 7r/60, much smaller than the 27r/m 
ambiguity in defining rotational phase from the wave- 
forms. 

Note that the early part of the (2, 1) waveforms con- 
tain, by far, the longest wavelength radiation present, 
suggesting a greater potential for problems caused by 
extracting the waveforms too close to the source, not 
yet in the wave zone. For the relatively deviant (2, 1) 
and (3, 2) modes, we performed Richardson extrapolation 
with respect to extraction radius i?oxt, using the values 
extracted at 45M and 90M and assuming an i?gxt error 
for each mode. These Richardson-extrapolated phases 
are also shown in the right panel of Fig. [T^ subtracted 
from the (2,2) rotational phase. Richardson extrapola- 
tion evidently reduces the early phase deviations in these 
modes considerably. 

The rotational phases calculated from the £ = m modes 
agree to within 0.025 radians during the inspiral for sev- 
eral hundred M before the peak E, effectively identical 
within the uncertainties of the numerical approach. This 
phase agreement is consistent with expectations based on 
the PN analysis, for which all (£,m) phases should agree 
up to 2PN order. These rotational phases also agree to 
within about 0.2 radians with the coordinate-dependent 
rotational phases measured from the puncture tracks, af- 
ter shifting the puncture track phase by an overall factor 
of 7r/2. Heuristically, we can think of each {£,m) ra- 
diation waveform mode as having been generated by the 
rotation of its own implicit source component. The phase 
agreement would then be interpreted as indicating that 
these implicit source components remain aligned through 
the inspiral. This is to be expected for a system which 
can be effectively described as an orbiting pair of point 
particles. 

It is, perhaps, more remarkable that a very tight agree- 
ment among the £ = m mode persists throughout the 
merger and even into the ringdown, remaining within 
about 1 radian until ~ 50A/ after the merger, when the 
amplitude has already diminished significantly. Accord- 
ing to our interpretation, this phase agreement suggests 
that a significant portion of the implicit rotating radia- 
tion source maintains some structural integrity through- 
out the coalescence. That is, the implicit source we 



have considered appears to exhibit considerable "rigid- 
ity" through merger. This is only possible because of the 
close relationship among the fundamental £ = m quasi- 
normal ringdown frequencies [see Sec. IV A| . mimicking 
the harmonic frequency relationship that holds during 
the inspiral. For the £ ^ m modes this quasinormal 
frequency relation docs not hold and the phases must 
separate in the merger. In terms of our implicit source 
picture, these £ ^ m components of the source seem to 
shear away from the main source structure to rotate at a 
faster rate. 

Heuristically, the puncture motion is strongly tied to 
the rotation of the implicit source for most of the evo- 
lution. During this period, it is natural to think of the 
implicit source as an inspiralling pair of pointlike objects 
moving on timelike world lines. At late times the orienta- 
tion phase angle of the puncture track disassociates from 
the waveform rotational phase. The punctures veer away 
from the implicit source at a late times as they fall into 
the final black hole. At this point, though we can con- 
tinue to consider an implicit rotating radiation source, it 
no longer makes sense to think of that source as a pair 
of pointlike objects. 

Having compared the phases of different multipolar 
modes for the 4:1 case, we now consider how the phase 
evolution depends on mass-ratio. There are various rea- 
sonable approaches to comparing the phases among sim- 
ulations of the different mass-ratio cases. Having estab- 
lished above the rotational phase consistency for the dif- 
ferent {£, m) modes, we will compare only the dominant 
(2, 2) phases. An obvious approach is to compare phases 
directly against time, scaled by the total PN mass M. In 
the early-time well-separated limit, however, the leading- 
order PN analysis indicates that phases for different mass 
ratios should evolve at similar rates when time is scaled 
by the chirp mass 

Fig. [T3I shows the rotational phase computed from the 
(2,2) mode of strain rate for different mass ratios. In the 
left panel, wc align the rotational phases at an early time 
and scale time by the chirp mass A4. For this plot, we 
shift the rotational phases in (chirp) time so that at t = 0, 
the chirp frequency, which is rotational frequency multi- 
plied by chirp mass, is 0.033 and the rotational phase 
is 0. Following this approach, we would expect good 
phase agreement at times sufficiently early that only the 
leading-order PN effects are significant. However, for 
the late portion of coalescence that we have simulated, 
we find that the different mass ratios remain roughly in 
phase for several hundred Ai before and after t = 0, 
peeling away in order at late times. 

In the right panel we compare phases in a manner com- 
mon for numerical relativity waveform comparisons, we 
shift the curves in time so that each peak energy flux 
occurs at t = 0, and we rotate the phases so that the 
phases are at this time. We scale the time by M. For 
the equal-mass case M = 0.435275, and for the other 
mass ratios it is smaller, so all of the curves in the left 
panel are stretched by at least a factor of 1/A4 w 2.3 in 
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FIG. 12: Rotational phase ^im calculated from puncture tracks and different multipolar strain-rate modes extracted at 45 M 
for a 4:1 mass ratio. The rotational phase is computed from each mode using (|8} and from the angle of the vector connecting 
the two punctures. (The puncture tracks remain in the x-y plane.) The phases are aligned so that the peak of E occurs at 
t = for all phases. The left panel compares the different calculations of ^tm, while the right panel compares the differences 
between the phases and the rotational phase from the (2,2) mode . In the right panel, we add Richardson extrapolations of the 
(2,1) and (3,2) modes based on the 45 M and 90 M extractions of these modes and assuming an error. The Richardson 
extrapolations are distinguished by asterisks in the figure. The differences are smoothed to better show the trends. 
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FIG. 13: Rotational phase $22 for different mass ratios, as computed from /i22. In the left panel, we time shifted the curves 
such that the chirp frequencies f222A^ ~ ^22M are equal at t = 0. In the right panel, we aligned the curves such that the /122 
peaks occur at t = 0, and the phases are at this time. 



time relative to the curves in the right panel. In the M- 
scaled right panel, the different mass ratios again remain 
approximately in phase for several hundred M before and 
after i = 0. At sufRciently late times, and particularly for 
small 77, we might expect this manner of consistent phas- 
ing as the evolution of the system eventually [after the 
innermost stable circular orbit] may become dominated 
by the course of unstable geodesic trajectories around 
the larger black hole [or an effective black hole in the 
effective-one-body (EOB) framework]. In that case the 
frequency evolves independently of the more strongly -q- 
dcpcndent rate of energy or angular momentum loss. 



V. DETAILED LATE-TIME ANALYSIS 

In Sec. lIVI we have presented general information about 
the phasing and amplitudes of the radiation components. 
Our analysis has stayed close to the standard numerical 
relativity waveform analysis, though we have emphasized 
an interpretation in terms of an implicit rotating source 
model. In this section we go beyond the standard wave- 
form presentation, exploring the radiation with the hope 
of developing a deeper understanding of the simple char- 
acteristics of the radiation as described above. Those fea- 
tures and our heuristic interpretation suggest a new ap- 
proach to examining the structure of the late-time phas- 
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FIG. 14: Rotational frequency for several {£, m) modes of the 
4:1 mass ratio simulation. These are time derivatives of the 
curves in Fig. 1121 For £ = m the rotational frequencies remain 
similar from inspiral through ringdown. These are similar to 
the puncture track frequency until t ~ — 40M, while the (2, 1) 
mode begins to deviate somewhat eailier. 

ing and apparent relationships between frequency and 
amplitude evolution. 

In Sec. IV A( we examine the phasing again, seeking 
a quantitative understanding of the late-time evolution 
of the polarization frequency. We introduce a practi- 
cal model that captures the merger-ringdown transition 
without the need for multiple quasinormal mode over- 
tones. We investigate the implications of this model in 
relating the frequency and amplitude close to merger in 
Sec. IVBI 

This section is more technical than Sec. IIVI with some 
subtle discussion of late-time radiation characteristics. 
For readers who may wish to jump ahead to Sec. IVIl 
we note two results that we will carry forward: (1) a 
simple quantification of the peak chirp rate lo for £ = m 
modes, and (2) the idea that dJ/duj becomes approxi- 
mately constant at late times, which may serve as sum- 
mary of relationships between frequency and amplitude 
near the radiation peak. 



A. Waveform frequency evolution 

First we will study the phasing in the merger and ring- 
down in more quantitative detail by comparing the po- 
larization frequency evolution for each mode with a sim- 
ple empirical model. Based on our heuristic model, we 
interpret the polarization frequency as corresponding to 
the rotational frequency {£, m) component of the implicit 
source. 

Fig. [Ml shows the evolution of the rotational frequency 
for several (^, to) modes of the 4:1 mass-ratio case. Sim- 
ilar frequency evolution curves are a common feature in 
papers on numerical relativity waveforms [see, for exam- 



ple, Refs. [l,i,[ili]. These curves are time derivatives of 
the phase evolution curves shown in Fig. 1121 zoomed in 
on the late-time behavior, near the elbows in the phase 
curves. The striking similarity in phasing for the various 
£ = m modes implies similar frequency evolution, which 
has been noted in previous studies (36| . 

At late times, this similarity in frequency is made 
possible because of a special approximate relationship 
among the fundamental £ ~ m quasinormal ringing fre- 
quencies, that they are nearly equal after dividing by 
the azimuthal mode number to to get what we call the 
rotational frequency. This has been considered in [37| . 
which pointed to a connection between the quasinormal 
ringing frequencies and the frequencies of stable null or- 
bits of a black hole at the "light ring" with frequency 
r^LR = l/{a + M(r+/M)(3/2))^ ^he association extends 
to charged Kerr-Newman black holes and has been com- 
pared with recent preci se q uasinormal ringing frequency 
calculations in [l3| and [38| . Conceptually, this allows us 
to think of the rotational frequency of the £ ~ m modes 
at late times as corresponding to the rotational rate of 
gravitational perturbations orbiting at the "light ring" . 
This suggests a heuristic description of our implicit ro- 
tating source at late times as a gravitational distortion 
of the forming final black hole which predominantly re- 
volves around the black hole on null orbits at the light 
ring. 

Returning to Fig. [TH we note that the (3, 2) mode is 
different from all the others, showing two spikes, near 
t — 9M and 37M. Comparisons of waveforms extracted 
at different radii, and from simulations of different reso- 
lutions, suggest some sensitivity to extraction radius, but 
do not suggest that the features will vanish in more ac- 
curate simulations or with more distant wave extraction. 
These anomalies may be related to the unusual shape 
in the amplitude peaks noted in Fig. [9] above. We will 
discuss this mode's behavior further in Sec. IVIII 

For all other modes the frequency evolution follows a 
simple smoothly evolving curve, qualitatively similar in 
each {£,m) and mass-ratio case [see Fig. [TSl below] . In 
particular, we note that, except for small noise contri- 
butions, each curve shows that the frequency increases 
monotonically, ultimately saturating at a frequency set 
by the fundamental quasinormal ringdown mode. This 
monotonic frequency development is a universal charac- 
teristic of the radiation from inspiral, through merger, 
and up to ringdown. In the PN analysis of quasicircular 
inspiral, this characteristic makes it possible to describe 
the changing structure of the hardening binary as a func- 
tion of frequency instead of the more coordinate-specific 
separation. This allows us, for instance, to write the 
waveform amplitude as a function of frequency. 

Subsequently, we will assume monotonic frequency de- 
velopment throughout the coalescence process. This 
principle underlies our empirical curve fitting of the fre- 
quency evolution, allowing more quantitative analysis of 
the late-time phasing evolution. In Sec. IVBl we will fur- 
ther apply this idea as we study relationships between 
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late-time frequency and amplitude evolution. 

To produce an empirical curve for describing the late- 
time frequency evolution, we assume that each {£, m) 
mode has a monotonically increasing polarization fre- 
quency, which approaches the fundamental ringdown fre- 
quency (jjQNM at late times. The general expectation that 
the frequency decays exponentially toward the ringdown 
frequency suggests that we model frequency evolution 
based on the hyperbolic tangent function. 

Specifically, we will compare frequency evolution of the 
strain-rate waveforms with a model of the form f2(t) — 
g{t), where 

(9) 

This provides a curve that first grows exponentially, with 
e-folding time 6/(2k), from some initial frequency fii, 
then decays exponentially, with e-folding time &/2, to 
the final frequency f2f. The presence of the exponent 
K allows the early exponential growth rate to differ from 
that at late times. The early part of this model is a coarse 
approximation to the growth in frequency near the end 
of the inspiral. This approximation must, therefore, fail 
to fit the data if we look back sufficiently long before 
the time of peak radiation, but, as we will show, it pro- 
vides a fair approximation of the approach to the peak. 
The rate at which the frequency grows, the chirp-rate, 
increases to some maximum, then decreases to zero on 
approach to the final frequency fif . For a more meaning- 
ful paramctrization, wc set 

b ■ ( iV""" 

r!i = r!f--r!ofi + -] • (lo) 

With this choice, dg/dt will peak with value fig at time 
to. The model then depends on the five parameters k, 6, 
i7o, to and fif. 

As shown in Figs. [TS] and (THl we find that fits to the 
model g{t) provide an excellent approximation to the nu- 
merical data for the strain-rate rotational frequencies in 
most significant cases. For the (2, 2) modes of all mass ra- 
tios, and for all but one of the significant (£, m) modes for 
the representative 4:1 mass-ratio case, we find agreement 
within a few percent after t = — 20M, with the primary 
differences coming from apparent noise in the numerical 
simulations, within the uncertainties in the numerical re- 
sults. 

Fig. [15] shows the comparison for the (2,2) mode ro- 
tational frequencies for several mass ratios. A glance at 
the curves shows that the unequal-mass cases arc quite 
similar to the equal-mass frequency evolution, previously 
examined in Refs. 0, [s^ . The dominant difference for 
the unequal-mass cases is that the final frequency rjf de- 
creases with 77, consistent with the decrease in the spin 
of the final black hole produced. We expect Jlf to cor- 
respond to wqnm/'ti, where wqnm is the fundamental 
(n = 0) quasinormal ringing frequency for the specific 



(^, m) mode for a black hole with the appropriate spin. 
In the infinite-mass-ratio limit (77 — > 0), fif should corre- 
spond to half the Schwarzschild (2, 2) quasinormal mode 
frequency AftjQNM = 0.3737, indicated by the horizontal 
dashed line in the left panel of Fig. [151 

We show a few examples of polarization frequency 
curves, for subdominant modes in Fig. 1161 The (3, 3) 
mode is very similar to the (2, 2) modes shown in Fig. ll5[ 
as are the other £ = m modes (not shown). The (2, 1) 
mode is of similar shape, also well approximated by our 
fit. 

The quantitative fit results are summarized in Table 
im The error bars are based on statistical fit estimates, 
also incorporating the ranges of best fit results obtained 
by varying the fit range starting between t/M = —25 and 
— 15 and ending between t/M = 20 and 60. The final fre- 
quencies rif approached in the fit curves in Figs. [15] and 
I16l werc robustly determined by the fits within a fraction 
of a percent. The ilf frequencies from the fits in Fig. [T5I 
were applied in Tabic |V] to find final black- hole parame- 
ters consistent with those determined by conservation of 
energy and angular momentum. 

The peak in the chirp-rate VIq is a particularly sig- 
nificant quantity in determining the shape of curves of 
our general form g{t). As shown in Table [TTl our fits 
determine Hq to within a few percent in all cases. An 
interesting relationship is apparent among the fits for all 
the £ = 711 cases for all values of rj studied. In the last 
column of the table, we show the peak chirp-rate values 
scaled by the mass of the final black hole Mf and the 
final frequency f2f. In each case where £ = m we find 
iloMf/rtf ~ 0.02, consistent within the fit uncertainties. 
This scaling makes some sense, since the height of the rise 
in frequency through the final radiation burst is largely 
determined by Jlf, while the time scale over which this 
rise occurs seems to be similar when time is scaled by 
M[. In Sec. IVII wc will use this result to predict the 
phase evolution in an analytic waveform model. 

Our model for late-time frequency evolution ((§]) de- 
scribes exponential decay toward Q{ at an c-folding rate 
given by half our fitting parameter b. For all cases, the 
values of b are within about 30% of lOM. In some cases, 
the fits for b are rather sensitive to the initial starting 
time, varying by up to 20% or 30% in the (2,1) and 
(4,4) modes of the 4:1 case. At this coarse level, we note 
that the values for b are similar to the exponential decay 
rates for quasinormal ringdown mode amplitudes listed 
in Tabic IIVI We will consider this relationship further in 
Sec. |VB] below. 

The other parameters in our fit are k, relating to the 
shape of our fit curve at early times, and to, giving the 
time at which the frequency peak occurs. The parameter 
K is not very precisely determined; as we would expect, it 
depends sensitively on the starting time of the fit inter- 
val, since the early exponential frequency growth is only 
a coarse approximation of the expected behavior. The 
values for to show that the peaks in Q generally occur 
roughly AM before the total energy peaks at ipoak- As 
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FIG. 15: Strain-rate rotational frequency evolution in the merger for (2,2) modes of several mass ratios (dashed lines), with an 
analytic fit (solid lines) . The fit encodes a monotonically increasing frequency, which at late times decays exponentially toward 
the fundamental ringdown mode. The horizontal dashed line marks Milt = 0.18685, half the (2,2) mode QNM frequency for a 
nonspinning perturbed hole. The right panel shows the residuals, which are comparable to uncertainties in the numerical data 
after 20M before the peak. 




FIG. 16: Strain-rate frequency evolution in the 4:1 mass-ratio merger for several spherical harmonic modes, with an analytic 
fit as in Fig. 1151 The left panel shows the polarization frequency ujirn- The right panel shows the residuals of the unsealed 
rotational frequencies Qim = uJ im/m. 



TABLE II: Results of fitting Eq. All quantities are scaled by the final mass Mf = Mf,rad defined in Eq. ([CTjl . The 

parameter Q,{ can be related to the fundamental real QNM frequency, ojqnm, and has been used to extract two different 
estimates of the final dimensionless spin a of the hole - see Table |V] 



Mass ratio 


{l,m) 


K h/M{ mCloMf^ 


mQfM{ 


to/M 




1:1 


(2,2) 


0.7 ± 0.1 13.2 ± 0.3 0.0112 ± 0.0001 


0.528 ± 0.001 


-4.6 ± 0.2 


0.0210 ± 0.0002 


2:1 


(2,2) 


0.6 ± 0.1 12.3 ± 0.8 0.0104 ± 0.0002 


0.5023 ± 0.0006 


-3.6 ±0.2 


0.0207 ± 0.0004 


4:1 


(2,2) 
(2,1) 
(3,3) 
(4,4) 


0.5 ± 0.1 10.5 ± 0.2 0.0096 ± 0.0001 
1.0 ±0.5 12 ± 2 0.0120 ± 0.0005 
0.3 ±0.1 9.7 ± 0.2 0.0153 ±0.0001 
0.15 ±0.1 7.5 ± 0.7 0.0212 ± 0.0003 


0.4566 ± 0.0001 
0.421 ± 0.006 
0.730 ± 0.001 
0.991 ± 0.005 


-4.8 ± 0.3 
-3.9 ±0.8 
-3.3 ±0.1 

-1 ± 1 


0.0210 ± 0.0002 
0.0285 ± 0.0012 
0.0210 ± 0.0002 
0.0213 ±0.003 


6:1 


(2,2) 


0.5 ± 0.1 10.4 ± 0.5 0.0089 ± 0.0002 


0.4349 ± 0.0004 


-6.5 ±0.3 


0.0205 ± 0.004 
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was the case for the power peaks in Fig.[8l the chirp-rates 
of the different spherical harmonic modes peak at slightly 
different times. 

We have supplemented our general implicit rotating 
source picture with the additional idea that the rota- 
tional frequency for each mode grows monotonically, not 
only in the inspiral, but also through the merger and ring- 
down. Based on this expectation we have identified an 
analytic fit model for the late-time frequency evolution 
that precisely matches the data for all cases but the (3,2) 
mode. These fits provide a quantitative understanding of 
the late-time phasing yielding, in particular, a robust re- 
sult for the peak chirp rate CIq for all ^ = to modes. We 
will apply this information in Sec. I VII 

B. Late-time frequency and amplitude 
relationships. 

The last step in our waveform analysis is to consider 
relationships between the frequency evolution and the 
waveform amplitude. 

In the PN description of the quasicircular inspiral, the 
orbital frequency not only tells us the rotational rate, 
but can also serve as a label for describing the momen- 
tary state of the rotating object (in the inspiral case this 
means that we can reference the state of the system in 
terms of r{U)). The PN generalization of the quadrupole 
formula, describing radiation from the rotating system, 
then leads to an expression for amplitude as a function of 
frequency. Our description of the gravitational radiation 
suggests an implicit source rotating with monotonically 
increasing frequency as it continues to "harden" , as the 
system evolves smoothly into merger and ringdown. In 
this section we seek to further unify this picture of the full 
coalescence process, considering an analogue of the PN 
description of amplitude as a function of frequency that 
can describe the radiation in the merger and ringdown. 

In Sec. IIVBI we emphasized that the radiation power 
Eim provides essentially the same information as the 
strain-rate amplitude (jA6p . If the wave frequency is 
known, then the modal contribution to the total radia- 
tive angular momentum can similarly provide informa- 
tion about the gravitational wave amplitude. Eq. (|A9[) 
gives an expression for angular momentum flux in terms 
of wave amplitude and phase. The relation simplifies to 

Jim « -^iA(.^)' (11) 

up to 5PN order, indicating that the radiation carries 
maximal angular momentum J = E/D,, as we gener- 
ally expect for circularly polarized radiation. If we know 
how the rotational frequency evolves, we can derive the 
mode amplitudes Aim from the mode-by-mode relation- 
ship of either energy or angular momentum with fre- 
quency, Eimii^) or Jimi^)- 

To approach an understanding of the late-time rela- 
tionships between amplitude and frequency, we compare 



how the system's energy and angular momentum ap- 
proach their final values, with how the system's frequency 
approaches its final value. In Fig. [17] we examine the 
radiative loss of energy dE and angular momentum dJ 
as the coalescing system approaches its final quiescent 
state, comparing these with the remaining difference of 
gravitational- wave frequency from its late-time limit f2f , 
as determined in Table HIl As well as our standard final 
rotational frequency = $ defined in ([5|), we also show 
the frequency: based on the strain, Qf^ = $/j as defined 
in ^ . We have also rescaled the energy and angular mo- 
mentum by a constant, selected so that the value matches 
that of f2/i at the time of peak radiation power. 

Fig. [17] indicates a general correspondence between 
how the angular momentum approaches its final state, 
and how the gravitational-wave frequency approaches its 
final state. Of the several cases of spherical harmonic 
modes and mass ratios that we have examined in this 
manner, we show two examples: the (2, 2) case for the 
equal-mass simulation, where the evolution of angular 
momentum and frequency correspond most closely (left 
panel), and the (2, 1) mode from the 4:1 mass ratio sim- 
ulation, with the weakest correspondence (aside from the 
nonconforming (3,2) case) (right panel). The correspon- 
dence is closest, holding to a better approximation over 
a longer period of time, in associating dJ with Clh- The 
association with energy is slightly weaker. 

Though these plots suffer significantly from small mod- 
ulations in frequency that we have not resolved numer- 
ically, the results suggest an approximate relationship 
between rotational frequency and angular momentum 
in particular. If the frequency evolution is otherwise 
known, then the late-time evolution of angular momen- 
tum for each mode could be approximately described by 
J{t) — Jf + £^{Q{t) — Q{), where ^ is a case-dependent con- 
stant that we will not attempt to specify generically. As 
^ ~ dJ/dD, we may refer to it as the dynamical moment 
of inertia of the implicit rotating source. This late-time 
expectation creates the possibility of extending our PN- 
based understanding of angular momentum flux into the 
late time waveforms, giving us the additional information 
we need for a full (approximate) description of J{fl). 

This relation between frequency, angular momentum, 
and amplitude provides a connection between the peaks 
in modal radiation power, shown in Fig. 1111 and the peaks 
in chirp rate given in Table [TTJ Assuming dJ/dft = ^, a 
constant, yields 

E w (12) 

=> E Pi ^ (nn + ii^^ (13) 

Since the ^l{t) curve is steeply increasing, we would ex- 
pect the peak in E to be near the peak in D, but slightly 
delayed. Expanding (|13p linearly about to, when f2 = 0, 
we find that E reaches its peak at a time 

tpeak ^ to - ■ ■ ■ > to, (14) 
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FIG. 17: Relationship between angular momentum and frequency evolution approaching ringdown for the (2,2) mode of the 1:1 
case (left panel) and the (2,1) mode of the 4:1 case (right panel). We show the differences dE and dj from the final mass and 
angular momentum as measured from the radiation, together with the difference from the final ringdown rotational frequency, 
measured from h {dfl) and from the strain h (dQh)- We see similar evolution for the dominant (2,2) mode of each mass ratio 
simulation, and for all modes other than (3,2). In the left panel the angular momentum and energy (in geometrical units) have 
been scaled up by factors of 2.5 and 10, respectively, to roughly match the frequency residuals at t = 0. In the right panel the 
rescaling factors are 760 for angular momentum and 2000 for energy. The agreement of the curves, with these rescalings, is 
consistent with dJ22 / dQ,22 ~ 0.40 from t = 20 through the peak radiation. The more approximate agreement in the right panel 
is consistent with dJ2i/dQ,2\ = 0.0013, through the peak radiation. 



where 57o and f2o arc the frequency and its third time 
derivative at to. These can be evaluated from ^ and 
the fit parameters k, b, CIq, and fif in Table [Til we find 
that the values for ipcak ~ all lie in the approximate 
range (2.16, 3.17), in rough agreement with the direct fit 
for io in Table [U 

Conceptually, the monotonic evolution of frequency 
might lead us to hypothesize a relationship between the 
structure of the implicit rotating source and its rotational 
frequency, such that the changes in the structure of this 
radiating rotator will be associated with finite changes 
in frequency. Knowing that the frequency growth must 
be limited by the quasinormal ringing frequency implies a 
peak in the chirp rate. This leads to the expectation that 
finite changes in energy and angular momentum are as- 
sociated with finite changes in frequency, so that dE/dfl 
and dJ/ c?f2 approach constants once the evolution in fre- 
quency slows. The peak in radiation power might then be 
viewed as a consequence of the peak in chirp-rate. More 
quantitatively we find that dJ/dil seems to be roughly 
constant even before the evolution in frequency slows 
down (i.e. before t = 0). This correspondence will be 
applied in the next section to provide a model for ampli- 
tude evolution through the peak, based on information 
about the frequency evolution. 

If we could postulate the constancy of dJ/dfl, we might 
also apply that assumption to "explain" some of what we 
have seen above. In Sec. IV Al we noted a rough agreement 
between the timescale b in our frequency fitting curve ^ 
and the quasinormal ringdown amplitude decay rates for 
the corresponding quasinormal modes. Following the dis- 



cussion above, this relationship could be derived, in the 
—* wqnm/™ limit, from the constancy of either energy 
or angular momentum losses with respect to change in 
frequency. For constant ^ = dJ/d^l, Eqs. (|A6|) . ([5]) and 
(dH) imply that 

^ 167r£: w 16T:^flfl. (15) 

In the Q ilf limit, the strain-rate amplitude decays 
at the rate predicted from black-hole perturbation theory. 
Aim ^o^mexp(-t/T£,„), where T£,„ is the e- folding 
rate for the amplitude decay for the fundamental {£, m) 
quasinormal mode. In this limit, our frequency evolution 
fit model reduces to 

n{t) f7f - (J7f - a)e-^(*-*")/^ (16) 

Applying these limiting expressions for amplitude and 
frequency in psp in the limit ^ Jlf yields 

(Ao)2e-^*/- « ^^M^(17f - f^Oe-'**-*"^/', (17) 
6 

where the left-hand-side derives from the amplitude, and 
the right-hand-side from frequency. Ignoring the con- 
stant coefficients, this implies that b = t. 

If we adventurously assume the constancy of ^ on ap- 
proaching the ringdown, and expand the amplitude in 
powers of e = exp(— 1/6) = exp(— t/r), the implied am- 
plitude frequency relation might also provide more infor- 
mation about the amplitude evolution. Since ^(t), and 
consequently the right-hand-side of (|17p . contains only 
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even powers of e, the next term in the expansion for am- 
phtude should be 0{e^). This suggestion motivated our 
expansion (jC4p applied in fitting the late-time amplitudes 
in Sec. |llll 



VI. VARIATIONS ON THE EOB MODEL 

An approach to modeling black-hole binary radiation 
known as the effective-one-body (EOB) model has been 
presented in the literature [H, M, M,\M,\M,M,\M- 
The late-time waveforms in these models are based on 
a now-common description of the merger process as an 
epoch of radiation from spiraling particlelike trajectories, 
followed, in a sudden transition, by black-hole ringdown 
dynamics with waveforms described by a superposition of 
quasinormal frequencies. Our waveform analysis provides 
a complementary description of black-hole binary merger 
radiation that can be applied in an alternative late-time 
waveform model. 

A recent promising approach along these lines is the so- 
called pseudo-4PN (p4PN) EOB model [l^. This model 
extends the 3PN-accurate EOB metric with a term of 
4PN order with a tunable multiplier A. The phasing 
obtained from this expansion is combined with leading- 
order FN strain amplitudes (the restricted approxima- 
tion) to obtain mode-by-mode waveforms valid for inspi- 
ral. For the merger and ringdown, the phasing and ampli- 
tude are derived simultaneously from a superposition of 
quasinormal modes. For each (£, m) angular mode, the 
fundamental ringdown modes and a few overtones are 
summed in proportions as required for continuity with 
the late end of the radiation from the inspiral phase. The 
value of the p4PN multiplier A is then chosen to match 
the pre- and post-merger waveform portions, optimizing 
the agreement with full- numeric waveforms. 

In this section, we show that it is possible to develop 
variations on the p4PN EOB model that usefully encode 
some of the waveform phase and amplitude relationships 
we have described above. A key difference with the new 
variant is our prescription for the transition from inspi- 
ral to merger-ringdown radiation. In contrast with [lO|, 
for each {i,m) mode we consider the entire wave train 
as that of a slowly varying instantaneously rigid rota- 
tor, consistent with the dominant "circular-polarization" 
waveform pattern encoded in the radiation. The phase 
evolution might be thought of as arising from the rotation 
rate of the corresponding {£, m) source structure, which 
is continued through ringdown by continuous matching 
to a function of the form ([9]). The wave amplitude will 
be derived directly from expectations for the energy or 
angular momentum content of the radiation. 

We present two specific models of this nature. Model 1 
is based on exactly the same EOB-based prescription for 
inspiral-plunge trajectories in [l^, while Model 2 shows 
the effect of a slight variation in the underlying EOB 
model. In both variants, as in p^ . we derive the wave- 
form phasing directly from the EOB trajectories (with 



A = 60 for the strength of the p4PN term) up until some 
matching time, which we take simply as the time at which 
the (2,2) wave frequency is half the frequency of the fun- 
damental (2,2) ringdown mode. After this point we use 
our fit model ([9]) to describe the subsequent phase evo- 
lution. 

Recall that this model depends on several parameters: 
K, 6, l^o, ^Oj E^nd f7f. The results of our analysis in Sec. 
IVBI guide us in producing a fully specified model for these 
parameters. We take J7f = ujQnM/'ni from the funda- 
mental ringdown frequency wqnm of the radiation. For 
the time-constant for frequency decay h we use the fun- 
damental quasinormal mode amplitude decay time con- 
stant. While this is not clearly implied by our fits in the 
last section, it will lead to the correct amplitude fall-off, 
as specified below. For the strongest £ = m modes, our 
fits indicate f2o = 0.021a;QNM/"^-A^f ■ Lacking any bet- 
ter model we simply increase this by a factor 4/3 when 
^ ^ TO, roughly consistent with the higher value of Oq 
found for the (2,1) mode. For these models we derive 
the quasinormal modes using the fit for the final black 
hole mass and spin described in p^ . The remaining pa- 
rameters to and K are chosen to provide continuity, up 
to the second time-derivative of phase, with the direct 
EOB-based phasing at the matching time. 

Our prescription for the wave amplitudes differs from 
the restricted amplitude description applied in [lol |. As 
we showed in Fig. 1111 we can improve on the restricted 
amplitude approximation by using the leading-order PN 
expressions for waveform (^, m) mode amplitudes only 
to fix the partitioning of radiation power into angular 
modes. We set the total power independently, from the 
full-order EOB model description of the radiation power; 
the resulting waveforms are then energetically consistent 
with the EOB description of the dynamics, and also show 
better agreement with the numerical results (note that 
in this model there is no radiation in the nonrotational 
TO = modes). After reaching the matching frequency, 
we continue the amplitude evolution based on the as- 
sumption, suggested in Sec. IV Bl that the amplitude is 
roughly constant through the radiation peak. 

In our new model, we set the late-time amplitude by 
asserting the approximate relationship (jlSp . written this 
time in terms of the polarization frequency uj = mU,: 

Am « 167r^£m (18) 

setting the value of for amplitude continuity at the 
match-frequency. With this model for the amplitude, 
the peak in the gravitational-wave amplitude is a direct 
consequence of the peak in the time-derivative of gravi- 
tational frequency, fixed by D,q in The exponential 
decay in amplitude also follows directly from the expo- 
nential approach of the wave frequency to wqnm • 

We compare the frequency and amplitude of the mod- 
eled waveform (2, 2) component with the corresponding 
numerical result for the 1:1 and 4:1 mass ratio cases in 
Fig. 1181 A similar comparison is shown for some of the 
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FIG. 18: Comparisons of new EOB-based models with numerical results for the (2,2) waveform frequency and amplitude for 
the 1:1 (left) and 4:1 (right) mass ratio cases. The models are variations on the p4PN BOB waveform model with a flux-based 
determination of the wave amplitudes and an alternative, based on Eq. [O] to quasinormal mode summing for continuing the 
waveforms through ringdown. The vertical bar indicates our matching frequency where we transition from direct BOB phasing 
and flux. Models 1 and 2 correspond to different versions of the radiative flux, which primarily affects the consequent wave 
amplitudes; see the text for more details. 




FIG. 19: Comparisons of new EOB-based models with numerical results for the 4:1 case, for the (3,3) (left) and (2,1) (right) 
waveform modes. The plots show frequency and amplitude, labeled as in Fig. 1181 



subdominant modes in the 4:1 case, in Fig. [121 The model 
we have so far described, based on the p4PN EOB tra- 
jectories is labeled Model 1 in the figures. The matching 
frequency is indicated as a vertical line in each plot. The 
frequency curves indicate very good phasing agreement 
for all cases except the (2, 1) mode of the 4:1 mass ratio 
case. In that case the sharp rise in frequency occurs a 
few M too late, several times worse than the agreement 
shown for the other modes. The difference is a conse- 
quence of the slightly higher frequency early on, as the 
numerical (2,1) frequency begins to grow already before 
the matching point as compared with the rotational fre- 
quency consistent with the numerical £ = m modes. 



The generally good amplitude agreement shown be- 
fore merger represents an improvement over the simple 
restricted amplitude model employed in [l3|. For the 
£ = m modes. Model 1 overestimates the amplitudes just 
before and consequently after matching. This indicates 
that the PN-based flux expression applied to the p4PN 
EOB model, based on [i^ and [i^l, overestimates the 
flux at high frequencies. This is perhaps not surprising, 
since the flux formula is formulated with a pole at finite 
frequency. Physically we expect the flux to decrease at 
late times, when the frequency approaches the quasinor- 
mal ringdown frequency. 

To correct for this, we show also results for a second 
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variation. Model 2. in which we have introduced a zero 
in the flux function at ringdown frequency. Following the 
notation of [i^ , we modify their Eq. (45) for the flux to 
read 

^Pm = -rV V — .fpN[v;v), (19) 

5 1-f/Wpolo 

and then respecify the coefficients ci — cj in their Eq. 
(50) to again provide consistency to 3.5PN order with 
the Taylor series expansion for the flux. Note that the 
flux now depends on £ and m via the quasinornial mode 
frequency wqnm and v = (MQ)^^^ ~ (A/w/m)^/^. The 
modified flux function anticipates that the radiation will 
cut off at the ringdown frequency. We find that the new 
EOB model provides a very good approximation to the 
original p4PN EOB phasing (with A = 60) if we choose 
A = 27 for the new version with modified fiux. The 
figures show that, for the i ^ m waveform modes, the 
amplitudes based on Model 2 with the modified fiux show 
better agreement with the numerical results leading up 
to the matching point, and also at the peak, while the 
frequency evolution is nearly identical. For the (2,1) case. 
Model 2 suffers the same problems as Model 1. 

We have applied the observations made in previous 
sections to successfully predict late-time waveforms. This 
provides an alternative description of the transition to 
ringdown, distinct from the widely-applied approach of 
summing quasinormal modes @, [111, SI] • In Ref . [l3| , we 
applied it to nonspinning mergers. 

This raises the question: how is it that each of these 
quite distinct approaches to approximately describing the 
late-time radiation can be simultaneously effective? We 
can explore this question by considering the polarization 
frequency evolution in a waveform constructed from a 
sum of quasinormal modes (suppressing £ and m labels) 



s{t) = ^yl„e-'""*+'^"(*\ 

n=0 



(20) 



with ipn{t) = </'n(0) + LUnt and where cun and (t„ = 1/t„ 
correspond to the n'^ quasinormal overtone mode. Next 
we restrict to the first two terms and to linear order in 
e(t) = exp[— ((Ti — cro)t], which vanishes at late times. 
This yields 

s{t) w ^oe-'""*+'"^«(*' (n-ee'(^i(*)-'^°(*») (21) 

~ ^pg-'^ot+e cos (Ai/>)+j[ipo+e sin (Ai^)] i'22) 

where Act = cti — ctq and A</j(t) = ipi{t) — ipo{t). Taking 
the derivative of the expression in square brackets gives 
the polarization frequency of s(t) 

ijJs (t) ^ ujQ — e [Act sin {A(p) — Aoj cos {A(p)] , (23) 

where Aw = lui —luq. Note that the expression in brackets 
is periodic with period 27r/Aw, where Aw is the differ- 
ence between the fundamental quasinormal ringing fre- 
quency and its first overtone [see [3] for a table of these 
overtones! . 



For the waveform modes we consider here MAw is 
quite small, generally 1 or 2%, which means that the pe- 
riod of the expression in brackets is > lOOM. We rewrite 
LUs, only keeping terms linear in Alu, as 



uJs{t) sa Wo — eACTsin (Aip) 



— — cot (Aip) 
Act ^ ^' 



(24) 



Considering the second term, we note that MAct is gen- 
erally just under 0.2 for the cases we've studied so that 
(Aw/Act) < 0.1. The other factor, cot(A(^), is not pre- 
dictable without knowledge of the initial conditions, but 
most values of Aip leave the second term somewhat less 
than one. 

If these conditions hold near the onset of ringdown, we 
might neglect the time dependence in the second term. 
The late-time frequency evolution is approximated by a 
exponential decay with a time constant of 1/Act. Work- 
ing under these assumptions, and comparing with the 
late-time frequency evolution model in this paper, yields 
the association 26 = 1/Act, where h is the fitting param- 
eter in our late time frequency evolution model. Again, 
looking at the quasinormal mode values, we comment 
that cti/cto ~ 3 so that 1/(2 Act) ~ & ~ tq, the approx- 
imate relationship noted in Sec IV A| which was also ap- 
proximately derived from the assumption that dJ / duj is 
constant at late times. Consistently, applying the rela- 
tion cti/cto ~ 3 in Eq. (|22p also give an expression for 
late-time amplitude consistent with Eq. (jC4|) . 



VII. DISCUSSION 

With hope of reaching out to a wide range of re- 
searchers interested in gravitational radiation from black- 
hole binary mergers, we have provided a descriptive walk- 
through of many of the general features of late-time wave- 
forms from generic mergers of nonspinning binary black 
hole systems, based on a series of numerical simulations 
covering systems from equal- mass up to mass ratio 6:1. In 
this basic waveform description we have examined wave- 
form phase and amplitude, comparing results among dif- 
ferent mass ratios, as well as among the different spin- 
weighted spherical harmonic (£, m) component modes. 

In our presentation, we have attempted to describe 
the radiation in the simplest physical terms, pointing out 
traits in the waveforms that are similar through the in- 
spiral, merger and ringdown stages. Throughout the co- 
alescence, we find simple waveforms in each (^, to) mode, 
each exhibiting strong circular polarization and mono- 
tonically increasing polarization frequency. 

In our amplitude comparisons, we find that the 
leading-order PN-prediction for energy-partitioning pro- 
vides a good estimate of the amplitude until late in the 
merger for £ = m modes. In astrophysical units, our fit 
for the peak-power in the dominant (2, 2) mode is 

E2Mti) ^ (1-607^2 + 19.7077"^) X lO^ergs-^ (25) 
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Scaling the amplitudes by the peak values yields a very 
similar shape through the peak for the (2, 2) mode ampli- 
tudes for all mass ratios we have studied. For a particular 
mass-ratio case, the peak-widths remain similar, though 
there is some variation in peak-time. 

For each mass ratio, the phase (and frequency) of the 
different {£, m) components are strongly related. While 
this should be expected for the early inspiral, where the 
waveform phase is directly connected to the orbital phase, 
we show that, for £ = to modes, the same relationship 
holds through the merger and into the inspiral. We com- 
pare the phasing among simulations with different mass 
ratios in two ways, with time scaled by chirp-mass, as is 
appropriate in the early-time limit, and then with time 
scaled by total system mass M . With the latter scaling, 
the dominant (2, 2) waveforms are similar in phase (and 
77-scaled amplitude) for the last ~ 200 A/. 

In the near-peak waveform comparisons the (3, 2) 
mode does not exhibit the same simple behavior as the 
other modes. This mode is easily subject to coupling 
with the much stronger (2, 2) mode. It has been seen 
in |3^ that the (3,2) mode demonstrates significant 
mode mixing with the dominant (2,2) mode - the QNM 
ringdown part of the (3,2) waveform contains the funda- 
mental frequencies of both the (3,2) and the (2,2) modes. 
We speculate that this mixing might be partly due to the 
use of coordinate extraction spheres that are systemat- 
ically warped from the areal-radius spheres appropriate 
for correct radiation extraction. A more refined choice 
of extraction spheres and perhaps a better tuned decom- 
position basis (e.g. spheroidal harmonics) would make 
it possible to represent the (3, 2) waveform content in a 
manner which is like that seen with the other modes. 
Thereby we expect that a similar simple physical repre- 
sentation of the waveform content can be extended to all 
TO 7^ modes. 

We suggest a simple conceptual interpretation that ap- 
plies through the full coalescence. We think of the radi- 
ation as being generated by an implicit rotating source^ 
with each {i, to) mode generated separately by the [i, to) 
moment of some implicit source (which we understand 
here only in the context of the radiation). The nearly 
fixed relationship among the (£, to) phase moments is 
interpreted to indicate that the implicit source main- 
tains some structural integrity throughout the coales- 
cence, without shearing among the various modal compo- 
nents. For the £ = m modes, this rigidity is maintained 
through the merger and into the ringdown, a relation- 
ship made possible by the approximate equality for each 
of the £ = m quasinormal modes cjqnm/™ ~ f^LR, where 
f^LR is the orbital frequency of unstable circular prograde 
graviton (or photon) orbits. 

The following physical picture may underlie these rela- 
tionships. For well-separated black holes, the fields that 
embody the implicit source object evidenced in the ra- 
diation may be tied directly to the pointlike centers of 
the orbiting black holes. The source rotation frequency 
is the orbital frequency of the timelike trajectories traced 



out by the black holes. As the binaries spiral together, 
the pair can continuously be viewed as a shrinking, dis- 
tributed dumbbell-like rotator. Eventually, most of this 
dumbbell shrinks inside the light-ring, which roughly co- 
incides with the potential barrier in the wave mechanics 
of gravitational perturbation theory. From inside, little 
radiation can escape to a distant observer, and the time- 
like motion of the black hole centers disconnects from the 
radiation. At late times, the effective radiation source be- 
comes a gravitational disturbance orbiting the forming 
black hole at the light ring. This is a seamless transi- 
tion, with nearly consistent rotational phasing among all 
£ = m modes throughout the process. For the £ ^ m 
modes the associated quasinormal-ringing dynamics are 
somewhat distinct, and the phasing and amplitude rela- 
tionships begin to peel away from the main £ — m trend 
through the merger process. 

For the late-time portions of the waveforms, includ- 
ing the approach to the peak and the ringdown, we have 
introduced a quantitative fitting model based on a mono- 
tonically increasing polarization frequency for each mode, 
which decays exponentially toward the expected funda- 
mental quasinormal ringdown frequency wqnm at late 
times. These fits provide an excellent match for the 
frequency evolution beginning ~ 20M before the peak, 
and allow precise estimates of wqnm , as well as the peak 
rate of change in frequency lvq. Scaling the latter quan- 
tity by the final black hole mass Mf and wqnm, we find 
tjo-A'jff/wQNM ~ 0.021 for all ^ = TO waveforms we have 
looked at, including the (2, 2) modes of each mass-ratio, 
and modes up to (4,4) for the 4:1 mass-ratio case. 

Conceptually, the monotonicity of the frequency evo- 
lution suggests that, as is the case for inspiralling sys- 
tems, the frequency can be taken to label the state of 
the adiabatically changing implicit rotating source that 
we interpret as the source of the radiation. If we sup- 
pose that changes in the source structure are tied to 
loss of angular momentum, then we would expect that fi- 
nite changes in frequency would be associated with finite 
angular momentum loss, so that dJ{u!)/doj has a finite, 
nonzero value even at late times. Since co approaches a 
nonzero constant at late times, we would likewise expect 
dJ/duj to approach a constant value. In Sec. IV Bi we 
show that the late-time evolution of J and lo are approx- 
imately consistent, mode- by- mode, with constant dJ/duj 
beginning about ^ 20Af before merger. 

Such a relation between frequency and angular mo- 
mentum also implies a connection between frequency and 
amplitude. The moment of peak amplitude is expected 
to be near the peak in coit). At very late times, constant 
dJ/dio implies a connection between the rates at which 
the frequency and amplitude approach their quasinormal 
late time state, namely that our fitting parameter 6 w r, 
where r is the damping time of the quasinormal ampli- 
tude decay [see Eq. [5] . 

The simple relationships between the waveform modes 
and simple dependence on mass-ratio make it possible to 
specify much of the late-time waveform information de- 
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veloped in our numerical simulations in terms of just a 
few quantities. This information can then be combined 
with information from the PN approximation about the 
inspiral trajectories to provide analytical models for full- 
coalescence waveforms. In Sec. IVII we have applied the 
PN-consistent EOB-based p4PN trajectory model pre- 
sented in [Toj together with assumptions asserting sev- 
eral of the approximate waveform features observed in 
Sec. IIVI As in [l^ , early waveform phasing is derived 
from the p4PN EOB trajectories, with waveform am- 
plitudes based on our PN-based power partitioning, to- 
gether with the PN flux model. As the waveform ap- 
proaches the anticipated peak, we match to a waveform 
phasing model based on our fit model, with parame- 
ters specified according to the approximate relationships 
identified in Sec. IIVI and with amplitudes derived from 
fiJ/(i(jj=constant. The fits show excellent phase agree- 
ment with the numerical simulation results for the most 
significant i = m modes. A variation in the PN flux 
model that enforces that the flux for each {£, m) mode 
vanish as — > wqnm gives better late-time amplitude 
agreement. 

Lastly, we observe that our description of the late- 
time phasing and amplitudes provides a picture comple- 
mentary to another approach applied in several previous 
studies d, [13, [m. El] , which successfully treat the late- 
time waveforms as a sum of quasinormal fundamental 
and overtone modes for each (£, m) waveform compo- 
nent. This is motivated by the expectation that wave- 
forms from generic initially compact distortions of the 
forming black hole will quickly reduce to a sum of these 
quasinormal harmonics [1^ . In [l0| , we have shown that 
in comparisons with some of the runs presented here, 
for mass ratios up to 2:1, this assumption can lead to 
a predictive waveform model with similar accuracy to 
our alternative model presented in Sec. IVII As a link 
between the two approaches, we have shown that at late 
times the combination of the fundamental and flrst QNM 
overtones for a particular mode may, under reasonable 
circumstances, mimic the amplitude decay properties of 
our model. 

This work suggests several directions for further study. 
In the immediate future, we plan to assess the fldelity 
of the available nonspinning waveforms and models, and 
the impact of mass ratio on the overall detectability of 
the merger signal. We also plan to apply our implicit- 
rotating-source description as a baseline in analyzing fu- 
ture higher-precision numerical simulations. This might 
provide insight into understanding finer features of the 
merger physics, some of which could violate our simpli- 
fied description. Further understanding the anomalous 
(3,2) mode waveforms will be a first step in this direc- 
tion. We must also investigate whether our description of 
the merger radiation applies also to spinning black holes. 
It is plausible that even precessing systems might be an- 
alyzed in this way usin g a spherical harmonic basis that 
tracks the orbital axis [50j. This may make it possible 
to extend the analytic EOB-based waveform model pre- 



sented here to include spin effects. Including spins in 
such analytic models will be necessary for observational 
data analysis applications. 
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APPENDIX A: RADIATION EXTRACTION 

We extracted radiation using the outgoing Weyl scalar 
'04, defined as in [5l|, calculated with a symmetric tetrad. 
This is related to the complex gravitational-wave-strain 
h{t, x) via 

V'4(t,^) = -h+{t,x) +ihx(t,x). (Al) 

"04, the strain h, and its time-derivative h (which we 
call the strain rate) are all functions of time t, extrac- 
tion radius i?cxt, and polar angles 6, <j). As is customary 
in numerical relativity, we decompose the radiation into 
spin-(-2)-weighted spherical harmonic components: 

oo t 

= Cer,r{t,R)^2Yr{e,^). (A2) 

(=2 m=~e 

oo e 

h = him{t,R)^2Yr{9,cP), (A3) 

1=2 m=-e 

oo e 

h = he^it,R)^2Yr{0,(p). (A4) 

e=2 m=-e 

With this, two time integrations of our measured quan- 
tity, -04, yields the more familiar gravitational- wave 
strain h. 

To extract the radiation from a simulation, we define a 
series of coordinate spheres of different radii i?Gxt ; here we 
use extraction spheres having radii between i?oxt = 40M 
and i?oxt = lOOAf . We extracted the radiation in modes 
by integrating ip4 against different -2Yf™'{6, (f>) over these 
coordinate spheres, using fourth-order interpolation onto 
each sphere followed by Newton-Cotes angular integra- 
tion. 

The gravitational waves produced by the binary carry 
both energy and angular momentum. Overall, the rate 
of energy emission is given by an angular integral of the 



22 



squared strain rate over a coordinate sphere [see Eq. 
(5.1) of dl]: 



= lim (/) dn\h\ 

dt R^oo 167r ^ ' ' 



j.|2 



(A5) 



Then using Eq. (|A3p . wc can express the total energy 
flux (IA5D as a sum over modes: 



1 



= lim '2 



'toll 



£7n 



167r 



(A6) 
(A7) 



where we have used the strain-rate decomposition (O and 
taken the limit i? ^ oo to go from (|A6p to (|A7p . 

Similarly, the rate of radiation of the z-component 
of angular momentum can be expressed as a sum over 
modes 1531: 



Irn 



~dt 



em 



lim 



RHui{h,^h*J. (A8) 



Substituting the definitions ^ and ^ into expression 
()A8P for the angular momentum, and taking the limit 
R-^ CO yields 



Jim = ^A(;„iHe,n cos(m($f,„ - (A9) 



APPENDIX B: CONVERGENCE 

We carried out three runs of the 4:1 mass ratio model 
at different resolutions to study the convergence proper- 
ties of our simulations. For these cases, the mesh spacing 
of the finest grids (the ones including the smaller punc- 
ture) was taken to be /i/ — 3M/160 (low resolution), 
hf = M/64 (medium resolution), and hf = 3M/224 
(high resolution). To facilitate comparisons among these 
cases, the overall grid structure of the runs was kept the 
same. In this Appendix, we discuss the convergence prop- 
erties of the constraints and gravitational waveforms. 

In comparing our medium and high resolutions for the 
4:1 mass ratio case, the Hamiltonian constraint was found 
to be manifestly fourth-order-convergent in the dynam- 
ical strong-field, where the black holes move, evidently 
dominated by the expected fourth-order error from re- 
finement interfaces. The convergence falls off to an ap- 
parent rate closer to first order in the coarsest regions. 
This seems to result from stronger dissipation of high 
frequency noise. Though less noise is generated in the 
higher resolution simulations, a greater portion of it sur- 
vives propagation into the distant coarser regions. As 
shown below, this does not appear to affect the wave- 
forms, which are well-resolved in the wave zone. The mo- 
mentum constraint appeared to be at least second-order- 
convergent in the dynamical strong-field region, but also 
fell off to an apparent rate closer to first order in the 
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FIG. 20: Convergence of strain-rate amplitude for (2,2) (up- 
per panel) and (4,4) (lower panel) modes. The medium-high 
differences are scaled for both fourth- and fifth-order conver- 
gence. Three-level convergence is lost around 150A'f before 
merger (vertical dashed line). 



coarsest regions - from the wave-extraction region out- 
wards. 

In our waveform mode analysis, we have typically time- 
and phase-shifted the data so that the amplitude peak 
fell at time i = In Fig. [3l the real part of the 
(2,2) and (4,4) strain-rate harmonics is shown for all three 
resolutions. The agreement between the resolutions is 
then seen to be excellent. In Figs. |4] and [5l we showed 
the errors to be expected from the maquillaged strain- 
rate waveform data. 

However, the presence of considerable eccentricity in 
the binary makes it difficult to compare such time-shifted 
waveforms between resolutions and establish an unam- 
biguous order of convergence. In In Figs. [^D] and [^T] we 
show the amplitude and phase errors, respectively, for 
our three resolutions without time-shifting - that is, we 
plot the data starting from the initial time in each case, 
scaling the medium - high differences for fourth and fifth- 
order convergence. It is clear that we observe conver- 
gence between fourth and fifth order throughout the evo- 
lution until around 150A/ before merger, when the differ- 
ence in merger times among the runs becomes important. 
Additionally, it is possible that at this higher-frequency 
stage of merger, the lowest-resolution data is no longer 
in the convergence regime. We do note, however, that 
the rate of growth of amplitude difference between the 
medium and high-resolution runs is comparable here to 
that observed closer to the merger for time-shifted data 
[see Fig. 2]. From this we deduce that the medium and 
high-resolution runs are still in the convergence regime, 
with errors consistent with fourth-order convergence. 
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FIG. 21: Convergence of strain-rate phase for (2,2) (upper 
panel) and (4,4) (lower panel) modes. The medium-high dif- 
ferences are scaled for both fourth- and fifth-order conver- 
gence. Three-level convergence is lost around 150Af before 
merger (vertical dashed line). 



APPENDIX C: END STATES 

In this Appendix we discuss the state of the final black 
hole formed in our simulations. This can be measured by 
several independent means, which we compare, finding 
agreement to within 0.4% for Mf and 2.1% for a. Ta- 
ble mil we present "coarse" results for the simulations. 
The total radiated energy AE^ad is obtained by integrat- 
ing Eq. (jASP , and the total radiated angular momentum 
AJrad is obtained by summing over mode contributions 
(|A8p . The final mass and angular momentum of the post- 
merger black hole are calculated using 

Mf,rad = AfADM - Ai;,ad, (CI) 
Jf,iad ^ Jq — AJiad, (C2) 

where, by the symmetries of the current simulations, we 
only deal with the z component of angular momentum. 
For the 1:1 simulations, we only had the leading-order 
(2, ±2) radiation modes available, so our A Jj-ad estimate 
is significantly truncated, by as much as 11% (a con- 
servative error estimate based on the effect of similarly 
truncating mode contributions past £ — 2 for the 4:1 
case); we have marked this and derived values. We also 
quote the measured value of tpeak, the time at which E 
reaches its peak. All waveform and derived plots in this 
paper have been time-shifted by subtracting this time, 
as an approximate marker of the time of merger, unless 
otherwise indicated. 

Table IIIII also contains data about the common ap- 
parent horizon (CAH) of the merged binary: the time 
^CAH at which the CAH was first detected, the CAH's 
irreducible mass Min-^cAH, and its full (horizon) mass. 



obtained from the Christodoulou [53| formula: 

M'cAH = A/2,. + ^, (C3) 

where we use J = Jf,iad for the final hole's angular mo- 
mentum. 

We can use the data from Table IIIII to estimate the 
mass and dimensionless spin of the end-state Kerr black 
hole through different methods. Given the final mass es- 
timates Mf.i-ad and Mf^cAH, we can calculate a = Jf/Mf. 

Another means of characterizing the final black hole 
comes from studying the characteristics of the radiation 
after the peak. This is expected to be a sum of the 
hole's quasinormal modes (QNMs). We can determine 
the imaginary part of the mode frequency by fitting the 
waveform amplitude A;_,„ from ([7]) post-peak to a de- 
caying exponential. The damping time in this fit corre- 
sponds to the imaginary frequency. Specifically, we fit to 
a functional form: 

A{t) = ^oe"*/^«''" (l - Ce"2*/^«™) . (C4) 

The presence of the additional damping term, 
parametrized by C, allows more freedom for non- 
linear decay early in the ringdown. Note that this factor 
of 3 between primary and secondary damping times 
approximately mirrors the difference in damping times 
one sees between the first two QNMs of a Kerr hole 
[111. The results of this fit are given in Table [IVl The 
value of tqnmj together with an estimate of the real 
QNM frequency wqnm, uniquely determines the mass 
and dimensionless spin of the Kerr hole. We combine 
our present fit for tqnm with the real frequency wqnm, 
as presented in Table HIl to obtain (MgNMiAQNivi)- 

If we already have a robust determination of the Kerr 
mass, then we only need one component of the QNM 
frequency to lock down a. In Table [III we use cjqnm to 
obtain aQNM,iG, assuming Mqnm = Aff^rad- 

We have gathered all the estimates of Mf and/or a 
discussed above in Table fVl It is instructive to compare 
these values with recent predictions for the final mass 
and spin from analytic and numerical methods. Working 
with a subset of the data supplied here, and also invoking 
the test-mass limit, [l^ suggest the one-parameter fits: 

Mf/M = 1 - (1 - V879)f7 - (0.498 ± 0.027)77^(C5) 
a = %/l277 - (2.900 ±0.065)772. (C6) 

The numerical data from Table |V] also fits well the 
wholly numerical formula (3.17a) from [ll[, as well as 
the parametrized formula of js^ . 

We note that when discussing the mass Mf of the end- 
state Kerr hole in Sees. IV Al and later, we have consis- 
tently taken the radiation-derived value, Mf = Mf^rad, as 
results from Q indicate this has errors up to a factor of 3 
less than QNM-based mass. Since a far larger proportion 
of the initial angular momentum than mass is emitted 
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TABLE III: Results of radiation and apparent-horizon analysis of merger. The total energy flux E (|A5[) reaches its peak at 
time t-pcak/M. A_Erad and AJrad are the total energy and (2) angular momentum radiated during the simulation (the latter 
calculated as a sum over modes (jASj), resulting in flnal mass Mf^rad and angular momentum Jf,rad (For the 1:1 case, only 
the 1 = 2 modes were available, so the emitted AJrad will underestimate, and Jf,rad will overestimate, the physical results.). 
Mirr.cAH is the irreducible mass of the common apparent horizon, first detected at simulation time tcAH. Aff^cAH is the mass 
deduced from this and Jf .rad using the Christodoulou formula [s^j . 



Mass 


hf 




h 


eak/A/ A£,,d(xlO-2)/M Mt,,,i/M 


AJ,,d(xlO-2)/M2 Jf,,ad/M2 


tcAu/M Mi„,CAH/M Mf,CAH/M 


ratio 
























1:1 




60 




1303.7 


3.5934 


0.9548 


33.79 


0.6468 


1196 


2:1 


3M/160 


45 




627.7 


2.8306 


0.9606 


23.34 


0.5937 


581 


0.9063 


0.9637 


4:1 


3M/224 


45 




641.9 


1.4327 


0.9786 


12.78 


0.4615 


588 


0.9489 


0.9796 




M/64 


45 




652.6 


1.4262 


0.9786 


12.99 


0.4594 


599 


0.9489 


0.9793 




3M/160 


45 




677.9 


1.4102 


0.9789 


12.90 


0.4603 


610 


0.9492 


0.9797 


6:1 


Af/64 


45 




564.9 


0.9212 


0.9850 


07.83 


0.3666 


513 


0.9667 


0.9851 



TABLE IV: Fit results for tqnm from Eq. (fC4)l . tqnm is 
used for one of the determinations of final Kerr parameters in 
Table |V] 



Mass ratio 




Ao(xlO-2) tqnm/A/ 


1:1 


(2,2) 


31.3 ±0.1 11.68 ±0.01 


2:1 


(2,2) 


26.80 ±0.06 11.57 ±0.01 


4:1 


(2,2) 
(2,1) 
(3,3) 
(4,4) 


17.23 ±0.09 11.37 ±0.01 

4.1 ±0.1 11.42 ±0.11 

8.2 ±0.2 10.96 ±0.06 
2.85 ±0.05 10.88 ±0.06 


6:1 


(2,2) 


12.77 ±0.08 11.35 ±0.02 



during the inspiral and merger, there is greater uncer- 
tainty in the radiation-derived final spin than the final 
mass. For this reason, our preferred measure of the final 
dimensionless spin is the QNM-derived value SqNM.ro- 
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Mass ratio 


A/f,rad/M 




Mf,CAH/M 


A/f,QNM/A/ 


Sqnm 


flQNM.ro 


1:1 


0.9548 


0.7095 




0.951 ±0.002 


0.684 ±0.002 


0.689 ± 0.002 


2:1 


0.9606 


0.6434 


0.9637 


0.961 ± 0.002 


0.624 ±0.001 


0.624 ± 0.002 
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0.9786 


0.4819 


0.9796 


0.979 ± 0.001 


0.472 ± 0.002 


0.4710 ± 0.0004 


6:1 
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0.9851 


0.989 ± 0.002 
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0.377 ± 0.002 
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